Monte Carlo study of the spin-glass phase of the site-diluted dipolar Ising model 
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By tempered Monte Carlo simulations, we study site-diluted Ising systems of magnetic dipoles. 
All dipoles are randomly placed on a fraction x of all L 3 sites of a simple cubic lattice, and point 
along a given crystalline axis. For x c < x < 1, where x c ~ 0.65, we find an antiferromagnetic 
phase below a temperature which vanishes asi->i c from above. At lower values of x, we find an 
equilibrium spin-glass (SG) phase below a temperature given by fcsTs 9 ~ xed, where Ed is a nearest 
neighbor dipole-dipole interaction energy. We study (a) the relative mean square deviation A 2 , of 
\q\, where q is the SG overlap parameter, and (b) (,l/L, where £l is a correlation length. From their 
variation with temperature and system size, we determine T sg . In the SG phase, we find (i) the 
mean values (| q \) and (q 2 ) decrease algebraically with L as L increases, (ii) double peaked, but 
wide, distributions of q/(\ q ]) appear to be independent of L, and (iii) £l/L rises with L at constant 
T, but extrapolations to 1/L — > give finite values. All of this is consistent with quasi-long-range 
order in the SG phase. 
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I. INTRODUCTION 

The collective behavior of spin systems in which mag- 
netic dipole-dipole interactions dominate has become the 
subject of considerable attention. These systems are rare 
in nature, although some ferroelectricsji and magnetic 
crystals such as LiHoF4, an insulating magnetic salt, have 
been known for decades to be well described by models 
of magnetic dipoles^— Much of the renewed interest in 
systems of interacting dipoles comes from the experimen- 
tal realization of magnetic nanoparticlo^ arrays^ and of 
crystals of organometallic molecules^ In these systems, 
particles up to some thousand Bohr magnetons behave 
as single spins. When closely packed in crystalline ar- 
rangements, dipolar interactions between them may in- 
duce magnetic ordering!^ 

Anisotropy also plays an important role in ordering 
dipolar systems. The barrier energies E a that must be 
overcome by spins in order to reverse their direction are 
often somewhat larger than the relevant dipolar ener- 
gies Ed- Then, collective effects can be observed when 
thermal energies are not sufficiently large to completely 
freeze spins directions. Their main effect is then to force 
spins to point up or down along the easy magnetization 
axisfi^ Crystalline Ising dipolar systems (IDSs) are then 
reasonable models^ These systems are clearly frustrated, 
since two different dipoles give rise to magnetic fields at 
any given site that are not in general collincar. Not sur- 
prisingly, IDSs are very sensitive to their spatial arrange- 
ment. Early work by Luttinger and Tisza established 
which type of magnetic order arises at low temperature 
for IDSs in each of the cubic lattices ~ More recently, we 
have obtained similar results by much simpler methods J£ 
For instance, BCC and LiHoF4 like crystals are ferromag- 
netic ordered, but antiferromagnetic (AF) order obtains 
on simple cubic (SC) lattices. Competition between dif- 



ferent interactions brings about a more exotic magnetic 
order, known as "spin ice" ^ in diamond crystals. 

Whether disorder in IDSs, together with the geomet- 
ric frustration that comes with the dipolar interactions 
give rise to a thermodynamic spin-glass (SG) phase, is 
an interesting question^ Many experiments^ as well 
as numerical simulations^ have shown that assemblies 
of classical magnetic moments placed at random, such 
as in frozen ferrofluids and diluted ferroelectric materi- 
als, exhibit the time dependent behavior, such as non- 
exponential relaxation and aging^ that is expected from 
SGs. However, search for evidence for the existence of 
an equilibrium SG phase has been hampered by the ex- 
tremely slow relaxation that is inherent to these sys- 
tems. In recent papers, we have given numerical evi- 
dence that supports the existence of an equilibrium SG 
phase in IDSs with randomly oriented axes both in fully 
occupied^ and in partially occupied SC lattices^ 

Site dilution is a rather simple way to introduce dis- 
order in experimental realizations of IDS. Some early at- 
tempts to find a SG phase in Eu^Si'i-^S led to nega- 
tive results^ By far the most scrutinized system for the 
last two decades has been LiHo a; Yi_ 2: F4. In it, magnetic 
Ho 3+ ions are substituted, with little distortion, by non 
magnetic Y 3+ ions.— A strong uniaxial anisotropy forces 
all spins to point up or down along the same axis at low 
temperatures. This parallel-axis-dipolar (PAD) system 
orders fcrromagnctically a low temperature phase above 
x c ~ 0.25. Below x c , transitions from a paramagnetic 
to a SG phase have been reported^— , but the oppo- 
site conclusion, that no such transition takes place, has 
been reached in Ref. [24]- The issue is further obscured 
by quantum effects that may take place at x <C 1^ 

Theoretical results suggest that diluted PAD models 
undergo a SG transition at low concentrations. An ear- 
lier study of bond-diluted Ising systems with long-range 
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interactions (including the dipolar case) found that SG 
order may exist at low temperatures in the limit of weak 
concentration.— Mean field calculations for site-diluted 
PAD systems in FCC and BCC lattices predicted a SG 
phase for concentrations < x < x c where x c is the value 
above which ferromagnetic order ensues^ More recently, 
Edwards- Anderson^ type models with power-law decay- 
ing interactions ~ ^/ r ij have been studied] 29 ' 30 A ID 
Ising Spin Glass model has been found to have a nonzero 
temperature SG phase transition for a < 1.— A 3D Ising 
systems with RKKY interactions (that decay with l/rf^) 
have been predicted to lie in the same universality class as 
the 3D Ising Edwards-Anderson (EA) model with short 
range interactions.— 

Numerical methods have provided conflicting answers 
to the question of the existence of a SG phase in site 
diluted PAD models. Biltmo and Henelius 3 -^ have calcu- 
lated that the ferromagnetic phase of LiH02.Y1_.jF4 ex- 
tends down to x c ~ 0.24, but found no SG phase at low 
temperatures for x < x c Q This is in contradiction with 
another MC simulation for the same system that finds a 
SG transition for concentrations x = 0.065 and 0.125.— 
Numerical work has also been done on a PAD model on 
a SC lattice, using a Wang-Landau MC method^ 3 - No 
transition was found for x < 0.2. 

Here we also simulate a PAD model on a SC lattice. 
Our justification for working with a SC lattice is as fol- 
lows. Whereas such systems order AF in fully occupied 
SC lattices , 11 ' 12 instead of ferromagnetically, as in the 
LiHoY4 lattice, the physics of PAD systems is not ex- 
pected to depend on lattice structure for x <C 1. A con- 
tinuum should then lead to the same behavior. Further- 
more, rescaling distance r as r — > rj p 1 / 3 , where p is the 
spatial density of spins, is no different from redefining 
dipolar energies by Sd — > psd, since dipolar interactions 
decay as r~ 3 . Now, consider kBT sg /rid£d for any lattice 
structure, where fcs is Boltzmann's constant, T sg is the 
SG transition temperature, rid is the number of magnetic 
dipoles within a d 3 volume, and Ed is the smallest possi- 
ble dipolar energy two parallel dipoles that are a distance 
d apart can have. Clearly, kBT sg /rid£d must be indepen- 
dent of lattice structure for x -C 1. This enables us to 
compare results for SC and LiHoF4 lattices, or any other 
lattice, for x -C 1. Such a comparison is made in Table 
I. 

The main aim of this paper is to find, by means of 
MC simulations, whether an equilibrium SG phase ex- 
ists in site diluted systems of dipoles, which are placed 
at random on the sites of a SC lattice and point up 
or down along a chosen principal axis. Since in the 
limit of low concentrations details of the lattice are ex- 
pected to become irrelevant, our results have direct con- 
nection with the experimental and numerical work men- 
tioned above. In this regard, we follow along the lines 
of Ref. [H. But we aim to go further. It is our pur- 
pose to also find whether the SG phase of the PAD 
model behaves marginally, that is, it has quasi-long- 
range order (as the XY model 3 -- in 2D), or whether it 



TABLE I: Spin-glass transition temperature for PAD systems. 
NIL is entered where a transition has been concluded not to 
take place. For LiHo a; Yi_ :E F4, we let d = 5.175 A, hence the 
mean number of spins in volume d 3 is n d = 1.926a: (since unit 
cells of LiHoYF 4 are 5.175 x 5.175 x 10.75 A 3 large and have 
4 Ho ions each 3 -); furthermore, Ed = 0.214 K.— On simple 
cubic lattices, we let d = a, hence rid = x. X3 is the nonlinear 
susceptibility, and v is the critical exponent for the correlation 
length. 



On LiHoYF4 type lattices 



Ref. 


Method 


X 


rid^d 


kBT sg /n d e d 


V 


21 


X3 


0.167 


0.069 K 


1.9 




22 


X3 


0.045 


0.019 K 


2.3 




23 


X3 


0.167 


0.069 K 


3.1 




24 


X3 


0.165 


0.068 K 


NIL 




24 


X3 


0.045 


0.019 K 


NIL 




31 


MC 


0.06 


0.025 K 


NIL 




31 


MC 


0.12 


0.049 K 


NIL 




32 


MC 


0.125 


0.052 K 


1.8 


1.3 


32 


MC 


0.0625 


0.026 K 


1.6 


1.3 



On simple cubic lattices 



Ref. 


Method 


X 


kBT sg /n d E d 


V 


33 


MC 


0.045,0.12,0.20 


NIL 




here 


MC 


0.35 


1.0(1) 


0.95 


here 


MC 


0.20 


1.0(1) 


0.95 



has spatially-extended states^- as in the droplet^- and 
replica-symmetry-breaking 3 -- pictures of the SG phase. 

The plan of the paper is as follows. In Sec. |n]we de- 
fine the model, give details on how we apply the parallel 
tempered Monte Carlo (TMC) algorithm? 3 - 9 , in order to 
get equilibrium results. We also define the quantities we 
calculate, including the spin overlap^ q, and £l, often 
referred to as a "correlation length" In Sec. Mil we 
give results for the dipolar AF phase we obtain for x > x c , 
where x c — 0.65, as well as for its nature and boundary. 
In Sec. I_V| we give numerical results we have found for 
(i) q distributions and (ii) £l/L, within the following x 
and T ran ges, 0.2 < x < 0.65 and 0.6a; < T < l.5x. In 
Sec. IV Al we examine the evidence we have in favor of 
the existence of a paramagnetic to SG phase transition 
when x < x c , and find that the transition temperature 
is given by fcgT S9 ~ x£d, where £d is a nearest neighbor 
dipole-dipole interaction energy which is defined in Sec. 
ITT1 In order to study the nature of the SG phase, we ex- 
amine the following evidence in Sec. IV Bl (i) the mean 
values (| q\) and (q 2 ) decrease algebraically with L as L 
increases, (ii) double peaked, but wide, distributions of 
q/(\ q |) appear to be independent of L, and (iii) £_/L 
rises with L at constant T, but extrapolates to finite val- 
ues as 1/1/ — > 0. We provide a specific example of spatial 
correlation functions which decay algebraically with dis- 
tance but lead to £_/L curves that spread out with L 
(for finite values of L) as T decreases below T sg , in rough 
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agreement with our MC results for £l/L. All of this is 
consistent with quasi-long-range order in the SG phase. 
In Sec. IV CI we find the best pair of values for T sg and 
v, to have curves £l/L f° r various values of L collapse 
onto a single curve if plotted vs (T/T sg — 1)L 1 / 1 ' over the 
T > T sg range. The values given in Table I are obtained. 



II. MODEL, METHOD, AND MEASURED 
QUANTITIES 

A. Model 

We consider site-diluted systems of Ising magnetic 
dipoles on a SC lattice. All dipoles point along the z 
axis of the lattice. Each site is occupied with probability 
x. The Hamiltonian is given by, 

« = 5EV^ (i) 

ij 

where the sum is over all occupied sites i and j except 
i = j, oi = ±1 on any occupied site i, 

T ij =e a (a/r ij ) 3 (l-3zt j /r* j ), (2) 

Tij is the distance between i and j sites, z,j is the z 
component of r^, e a is an energy, and a is the SC lattice 
constant. In the following we give all temperatures and 
energies in terms of e a /&B and e a , respectively. Hence, 
kBT/n a Sa = T/x from here on. 

This model is clearly an Ising model with long-range 
interactions where bond strengths Tjj are determined by 
the dipole-dipole terms. Note that TJy signs are not dis- 
tributed at random, but depend only on the orientation 
of vectors ry on a SC lattice. This is to be contrasted 
with a random-axes dipolar model, (RAD)i£ in which 
Ising dipoles point along directions iii = (nf,a = 1,2,3) 
that are chosen at random by sorting two independent 
random numbers for each site, introducing randomness 
on bond strengths T°f . This is why PADs exhibit AF 
order at high concentration in contrast with RADs, that 
do not 



B. Method 

We use periodic boundary conditions (PBC). As is 
usual for PBC, think of a periodic arrangement of repli- 
cas that span all space beyond the system of interest. 
These replicas are exact copies of the Hamiltonian and 
of the spin configuration of the system of interest. De- 
tails of the PBC scheme we use can be found in Ref. 
[l2l . We let a spin on site i interact through dipolar fields 
with all spins within an L x L x L cube centered on it. 
No interactions with other spins are taken into account. 
This introduces an error which we show in Appendix I to 
vanish as L — > 00 , regardless of whether the system is in 



TABLE II: Parameters of the tempered MC simulations, x is 
the probability that any given site is occupied by a magnetic 
dipole; L is the linear lattice size; AT is the temperature step 
in the TMC runs; T and T n are the highest and the lowest 
temperatures, respectively; N r is the number of (quenched) 
disordered samples; a number to of MC sweeps are made be- 
fore any measurements are taken. The measuring time inter- 
val is [to, 2to] in every case. 





X 


= 0.20, Ai = 0.02, io = 0.8 




L 


i 


6 8 10 




T n 


0.06 


0.06 0.06 0.12 




N r 


8500 


tr iar\r\ -\ r\r\r\ ar\r\ 

3800 1000 800 




to 


5 x 10^ 


O X 1U O X 1U O X 1U 






X 


= O.oo, AJ = U.05, lo = 2.0 




L 


1 


6 8 10 


12 


T„ 


0.05 


0.05 0.05 0.275 


0.35 


N r 


9000 


5000 1100 380 


200 


to 


4 x 10^ 


4 X 1U 4 X 1U 4 X 1U 


4 x 10^ 




X 


= 0.50, AJ = 0.05, lo = 2.0 




L 


1 


6 8 10 




T n 


0.1 


0.05 0.05 0.35 




N r 


1000 


n E c\ E f\r\ * 1 < 1 / -1 

650 500 300 




to 


5 x 10 


O X 1U 4 X 1U 10 






X 


= 0.60, AJ = 0.1, io = 2.0 




L 


i 


6 8 10 




T„ 


0.10 


0.10 0.20 0.30 




N r 


1400 


500 800 300 




to 


4 x 10^ 


4 X 1U 4 X 1U 4 X 1U 






X 


= 0.65, AJ = 0.1, Jo = 3.0 




L 


4 


6 8 10 




T n 


0.10 


0.10 0.10 0.30 




N T 


1400 


900 1400 540 




to 


4 x 10^ 


4 X 1U 4 X 1U 4 X 1U 






X 


= 0.70, Ai = 0.1, io = 3.0 




L 


4 


6 8 10 




T. 


0.10 


10 10 30 




Nr 


750 


200 100 100 




to 


4 x 10 6 


4 x 10 6 4 x 10 6 10 6 






X 


= 0.75, AT = 0.1, T = 3.0 




L 


4 


6 8 10 




T n 


0.10 


0.10 0.10 0.10 




N r 


1000 


200 100 100 




to 


4 x 10 6 


4 x 10 6 2 x 10 6 10 6 






X 


= 0.80, AT = 0.1, T = 3.0 




L 


4 


6 8 10 




T n 


0.10 


0.10 0.10 0.10 




N r 


600 


200 220 100 




to 


4 x 10 6 


4 x 10 6 10 6 10 6 





the paramagnetic, AF or SG phase. There is, therefore, 
no effect on the thermodynamic limit of the system of 
interest here. (The result we obtain in Appendix I is not 
applicable to an inhomogeneous ferromagnetic- phase or 
critical region- that may obtain on other lattices.) 

In order to bypass energy barriers that can trap a 
system's state at low temperatures in the glassy phase 
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we have used the parallel tempered Monte Carlo (TMC) 
algorithm . 39 ' 43 We apply the TMC algorithm as follows. 
We run in parallel a set of n identical systems at equally 
spaced temperatures Tj, given by Tj = Tq — iAT where 
i = 0, • • • , n — 1 and AT > 0. By identical we mean here 
that all n systems have the same quenched distribution 
of empty sites, though each system starts from an inde- 
pendently chosen initial condition. We apply the TMC 
algorithm to any given system in two steps. In the first 
step, system i evolves independently for 8 MC sweeps un- 
der the standard single-spin-flip Metropolis algorithm^ 
(Owing to dipolar interactions, the MC sweep time scales 
as N 2 , where iV is the number of spins.) We update all 
dipolar fields throughout the system every time a spin 
flip is accepted. In the second step, we give system i 
a chance to exchange states with system i + 1 evolving 
at a lower temperature Tj — AT. We accept exchanges 
with probability P = 1 if 5E = E- t — E i+ \ < 0, and 
P = exp(-Afi6E) otherwise, where A/3 = 1/T m -1/T;. 
The cycle is complete when i has been swept from to 
?i — 2. Thus, we associate eight MC sweeps with each 
cycle. For the simulation to converge at low tempera- 
tures it is important to choose AT small enough to allow 
frequent state exchanges between systems. This will of- 
ten be fulfilled if A/3AE < I. The required condition, 
AT < Tj \J Nc, follows for AT where c is the specific 
heat per spin. Then, we obtain appropriate values for 
AT from inspection of plots of the specific heat vs T»i£ 
We find it helpful to have the highest temperature To at 
least twice as large as what we expect to be the transi- 
tion temperature between the paramagnetic and the or- 
dered phase for obtaining equilibrium results in the or- 
dered phase. 

In our simulations the n identical systems start from 
completely disordered spins configurations. We need 
equilibration times to of at least 4 x 10 6 MC sweeps for 
x < 0.7 for systems with a number dipoles N > 200 (see 
at the end of this sections for details on how we choose 
to)- Thermal averages come from averaging over the time 
range [to, 2to]. We further average over N r samples with 
different realizations of disorder. Values of the parame- 
ters for all TMC runs are given in Table I. 



the moments 



(4) 



for n = 1,2, where (...) stand for averages over time 
and over a number N r of system samples with different 
quenched disorder. Unless otherwise stated, time aver- 
ages are performed over a time range to < t < 2t , and 
to is chosen as specified below in order to ensure equi- 
librium. We make use of these moments to calculate the 
staggered susceptibility and the mean square deviation 
of | m | /mi, that is, 



A 



m 2 



(5) 



In order to spot SG behavior, we also calculate the 
Edwards- Anderson overlap parameter^ 



where 



of of 



(6) 



(7) 



trj 1 ' and crj 2 ' are the spins on site j of identical replicas 
(1) and (2) of the system of interest. As usual, identical 
replicas have the same Hamiltonian, and are at the same 
temperature, but are in uncorrelated states. Clearly, q is 
a measure of the spin configuration overlap between the 
two replicas. As we do for m, we calculate the probability 
distribution P q as well as the moments q\ = (\q\) and 
q-2 = (q 2 ), in analogy to Eq. (f4]). The SG susceptibility 
Xsg is given by Nq 2 . Finally, we also make use of the 
relative mean square deviation of q, A 2 = q 2 jq\ — 1. 

We need to make sure that equilibrium is reached be- 
fore we start taking measurements. To this end, we de- 
fine a time dependent spin overlap q, not between pairs 
of identical systems, but between spin configurations of 
the same system at two different times to and t\ = to + 1 
of the same TMC run, 



q(t ,t) = N- 1 J2°j( t o)(T j (to + t). 



(8) 



C. Measured quantities 

We next specify the quantities we calculate. Wc obtain 
the specific heat from the temperature derivative of the 
energy. For the staggered magnetization, we define, as 
befits a PAD model on a SC lattice 1 ^ 

m = N- 1 ^2a i (-l) x M+y( i ) (3) 

i 

where x(i) and y(i) are the space coordinates of site i. 
We calculate the probability distribution P m , as well as 



Let q~2(t>o,t) = ([q~(to,t)} 2 }. Suppose thermal equilib- 
rium is reached long before time to has elapsed. Then, 
Q2{to,t) —> <Z2 at some time t long before t = t - Plots 
of q2(t ,t) vs t, for 10"% < t < t , for t = 10 7 MC 
sweeps, are shown in Fig. [T] for x = 0.20 and various 
values of T. Plots of q 2 , obtained by averaging q 2 over 
time, not starting at t = to, as we do everywhere else in 
order to obtain equilibrium values, but starting at t = 0, 
from an initial random spin configuration, are also shown 
in Fig. [T] for comparison. Note that both quantities do 
become approximately equal when t > 10 5 MC sweeps. 
In order to obtain equilibrium results, we have always 
chosen sufficiently large values of to to make sure that 
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FIG. 1: (Color online) Semilog plots of faito, t) and qi vs time 
t (in MC sweeps) for systems of 8 x 8 x 8 spins at the values 
of T shown in the figure. Here, 52 comes from averages of 
q 2 over time, starting at t = from an initial random spin 
configuration. Here, to = 10 7 MC sweeps. A data point at 
time t stands for an average over a time interval \t, 1.2t\ and 
over 10 3 system samples. 



Qzitoit) —> q2 long before t = to. All values of to and N r 
are given in Table II. 

As has become customary in SG work,— ~— we calcu- 
late quantity 



1 



4sin^(fc/2) 



9(k) 



- 1 



where 



g(k) = N' 



i 



(9) 



(10) 



Tj is the position of site j, and k = (2tt/L, 0, 0). Re- 
call this system is anisotropic, interactions along the spin 
axes are twice as large as in a perpendicular direction. 
We have found this direction of k (perpendicular to all 
spin directions) to be more convenient to work with than 
the direction along the spin axes. 



Note that replacement of exp(ik.rj) by 1 
E l7 [k.(r 2 - r,)] 2 (^) 



el 



8sin 2 (fc/2)£ 



ik..Yj gives 



(11) 



This is right in the £l/L — > limit. The above equation 
clearly shows that £t is then (up to a multiplicative con- 
stant) the spatial correlation length (in the k direction) 
of {(j>Q(j) r ). Therefore, we can think of ^ 00l the L — > 00 
limit of £l , as the correlation length of a macroscopic sys- 
tem in the paramagnetic phase. In a condensed phase, 
on the other hand, condensate fluctuations generally take 
place over finite lengths £, but (,l/L — > 00 as L — > 00 if 
there is strong long-range order, that is, if {4>o4> r ) does 
not vanish as r — > 00. One would have to replace <p by 
<j> — (cj)) in Eq. © in order to relate ^ to £. Following 
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FIG. 2: (Color online) Phase diagram of the PAD model, o 
stand for the Neel temperature Taf, and ■ stand for the SG 
transition temperature T sg . o stand for maxima value of x 
for which m2 decreases as N increases for each of three fixed 
values of T. The full line for the phase boundary between 
the paramagnetic and AF phases is a fit to the data points, 
given by, Taf — 3.8(x — x c ) 0A , where x c — 0.65. The straight 
dashed-line is for T ag — xe a - In the inset, versus x for 
T = 0.4. o, □, o, and A, stand for L = 10,8,6, and 4 
respectively. 



current usage, we shall nevertheless refer to £l as "the 
correlation length" . 

In contrast with P q and its first moments, £l takes into 
account spatial variations of the EA overlap q and is yet 
another probe for detecting a SG transition.— _ — 



III. THE AF PHASE 

Our main results for the PAD model are summarized 
in the phase diagram exhibited in Fig. [2l A thermally 
driven second order transition takes place at the phase 
boundary between the paramagnetic and AF phases. The 
phase boundary meets the T = line at x ~ 0.65. We 
shall refer to the value of x at this point as x c . 

In this section we report the numerical evidence for the 
paramagnetic- AF transition^ Results having to do with 
the spin glass are given in the next section. 

The AF phase is defined by the staggered magneti- 
zation, as given in Eq. ©. We illustrate in Fig. [3k. 
how the staggered magnetization mi behaves with tem- 
perature for x = 0.8. This is in sharp contrast to the 
behavior of mi for small x, where an AF phase does not 
exist. Such behavior is exhibited in Fig. (5}d. Note that 
mi appears to decrease as N increases even at low T. We 
obtain similar results for the staggered magnetization for 
other values of x (shown in Fig. [2|) below x c . This is 
our first piece of evidence for the nonexistence of an AF 
phase below some x c and that x c ~ 0.6. We return to 
this point in the discussion of Fig. 2J 

Plots of the specific heat C vs T arc shown in the 
insets of Figs. [3k. and [3J)- Note the sharp variation of 
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FIG. 3: (Color online) (a) Staggered magnetization mi vs T 
for x = 0.8. Icons o, □, o, and A stand for L = 10, 8, 6 and 4 
respectively. Lines are only guides to the eye. Note mi grows 
with L at low temperature, consistently with an AF phase. 
In the inset, specific heat vs. T for the same values of x and 
of system sizes. The sharp variation C with respect to T near 
T — 1.5 is consistent with an AF phase transition thereon, (b) 
Same as in (a) but for x = 0.6. Note (i) mi decreases with 
L at all temperatures, consistently with the nonexistence of 
an AF phase, and (ii) a rounded specific heat, consistent with 
a SG transition. In all panels, error bars are smaller than 
symbol sizes. 
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FIG. 5: (Color online) (a) Plots of A'f n vs T, for x = 0.7. •, 
□ , o, and x are for L = 10, 8, 6 and 4, respectively. Lines are 
guides to the eye. The thick dashed-line is for the macroscopic 
paramagnetic limit n/2—1. (b) Same as in (a) but for x = 0.6. 
(c) Plots of vs T, for x = 0.7. Symbols are as in (a), (d) 
Same as in (c) but for x = 0.6. Error bars are shown only 
where they are larger than symbol sizes. 




disordered phase, above T = 1.2(1), for which A^m 2 = 
0(1), to a strong long-range order below T = 1.2(1), 
where m 2 = 0(1). Note that m 2 ~ l /N p at T = 1.2. 
From the definition of r] (see Sec. IVBI or Ref. HH), 3p = 
l + ?7 follows, which gives r\ = 0.05. We are however not 
too interested here in such details of the critical behavior 
on the T = T AF (x) line. In Fig. [4}j, m 2 vs N plots 
show faster than algebraic decay with N . This shows we 
are then beyond the bounds of the AF phase. We have 
followed this criterion as a first approach in establishing 
the boundary of the AF phase. Plots of m\ (instead of 
m 2 ) vs show the same qualitative behavior. 

We draw more quantitative results about the AF phase 



FIG. 4: (Color online) (a) Log-log plots of m,2 versus iV for 
x — 0.7 and the values of T shown. Continuous lines are 
guides to the eye, except for the straight line over the data 
points for T — 1.2, which is for 1/A 0,35 . A dashed line shows 
the slope one expects for a macroscopic paramagnet. (b) 
Same as in (a) but for x — 0.5. In all panels, error bars 
are smaller than symbol sizes. 



C vs T near T = 1.5, in Fig. [3^, as one expects from 
a paramagnetic- AF phase transition. Note also how, as 
one expects for a paramagnetic-SG transition, C varies 
smoothly for a smaller value of x, in Fig. [3Jd. 

For further information about the extent of the AF 
phase, we now examine how m varies with N for some 
values of x and of T. Compare the log-log plots of ?ti 2 
versus the number of dipolcs N on Figs. 0^, and UJd, 
respectively. The data points in Fig. 0Ji are consistent 
with a second order phase transition from a magnetically 




FIG. 6: (Color online) Semilog plots of q2 versus T for x = 
0.2, and L = 10, (o), L = 8 (□), L = 6 (o), and L — 4 (>). 
All error bars are smaller than symbol sizes. 
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boundary from the behavior of the relative uncertainty 
Af n . We first outline how we expect to behave as a 
function of T and x in the various magnetic phases. It 
clearly follows from its definition in Eq. ([SJ that A, n — > 
as N — > oo in the AF phase. It also follows immediately 
from the the law of large numbers that, in the paramag- 
netic phase, — s- 7r/2 — 1 as TV — s- oo. These two state- 
ments imply that curves of A^ vs T for various values 
of N cross at the phase boundary between the paramag- 
netic and AF phases. We make use of this fact to quanti- 
tatively determine the AF-paramagnet phase boundary. 
The same criterion can be applied to the AF-SG phase 
boundary. To see why this is so, note that, the plots 
shown in Fig. [4}) for x = 0.5 suggest m,2 — > N^ 1 as 
N — > oo, even at low temperatures, that is, well within 
the SG phase. Plots of A^ vs T are shown in Figs. [5^ 
and [5}} for x = 0.7 and 0.6, respectively. The signature 
of an AF phase below T ~ 1.2 clearly shows up in Fig. 

We have thus established all points of the AF phase 
boundary shown in Fig. [5] for x > 0.7. For the low tem- 
perature portion of the phase boundary (near x = 0.65) 
this procedure is not very effective. From Fig. [SJd, we 
infer that the AF boundary line must drop to a T = 
value at some x > 0.60. The three data points shown for 
x ~ 0.65 and T < 1 are obtained from plots such as the 
one shown in the inset of Fig. [2] for T = 0.4. 



IV. THE SG PHASE 

In this section, we report numerical results we draw 
from tempered MC calculations for q2, for distributions 
of q, and for £j> Because we expect, from the argument 
given in Sec. I, lattice independent behavior for x <C 1, 
we emphasize the results we have obtained for the two 
smallest values of x we have dealt with, x = 0.2 and 
x = 0.35 (that is, x ~ 0.3a; c and x ~ 0.54a; c ). 

A plot of q2 versus T is shown in Fig. [6] Note that q^ 
decreases as N increases, even at low temperatures. We 
have found similar behavior for other values of x satisfy- 



1.00 



ing x 



< 



Inspection of this figure raises the question 
of whether q2 vanishes as L — > oo . In order to advance in 
this direction, we do log-log plots of q2 vs N, which we 
show in Figs. [Th., [7b, and [7b, for the values of x shown 
therein. The data points in these three figures seem con- 
sistent with, q 2 ~ N~ p for T/x < 1, where 3p = 1 + rj, as 
follows from the definition of r\ in Sec. EH (see also Ref. 
l4rf ). x 2 values for q2 ~ N~ p fits to sets of data points, for 
T/x < 1 (for which they are appropriate) as well as for 
T/x > 1 (for which they are not appropriate), are given 
in Table III. Plots of qi vs N show the same qualitative 
behavior. All of this is in accordance with quasi-long- 
range order. We return to this point below and in Sec. 

Reading off values of p from plots shown in Figs. [TJi, 
[7b>, and [7b, we obtain rj for x < 0.5 and various values of 
T. The relation 7] — — 1 + a x (T/x) 2 fits the data rather 
well for all T/x < 1, if we let a x = 0.76, 0.98, 1.18 
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FIG. 7: (Color online) (a) Plots of q2 versus the number of 
dipoles N for x = 0.5. o, □, o, >, A, V, <, and ■ stand for T = 
0.1,0.2,0.3,0.4,0.5,0.6,0.7, and 0.8 respectively. Lines are 
guides to the eyes, (b) Same as in (a) but for x = 0.35. o, □, o, 
>, A, V, and <, stand for T = 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, and 
0.6 respectively, (c) Same as in (a) but for x = 0.2. o, □, o, 
>, A, V, «, and ■ stand for T = 0.12, 0.14, 0.16, 0.22, 0.26, 
and 0.30. For all data, we have checked that, within errors, 
92 = 92- Clearly, data point sets for larger temperatures 
deviate from the straight dashed-lines shown (implying faster 
than a power of 1/L decay) while sets for lower temperatures 
do not. Error bars are shown only where they are larger than 
the icon sizes. For each set of points with given x and T 
values, x 2 values for straight line fits, as well as the largest 
error, are given in Table III. 



for x = 0.2, 0.35, 0.5, respectively. In order to be able 
to conclude that r](T sg ) varies with x, we would need to 
know T sg within an error of 10%. Unfortunately, we find 
below (in Sec. IV A[) an error in T sg which is not much 
smaller than 10%. 

For higher values of T/x, q2 vs TV curves downwards, 
as expected for the paramagnetic phase. Approximate 
values of T sg can thus be obtained from such plots, but 
more accurate methods are given below. It is reassuring 
to see in Figs. [7K, [7t> and [7fc, the values of q~2 we have 
obtained agree, within errors, with the values for q 2 . 

We next give distributions of q we have found. We 
make use of a normalized distribution P q (q r ), where 
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TABLE III: Xr values for two-parameter 52 = c/N p fits to 
sets of data points for 52 vs T displayed in Figs. [T^i-a As 
usual, we define Xr = X 2 /df, where df is the number of data 
points in each set minus the number of fitting parameters (2, 
here). The largest errors A52 of 52 from all data points for 
each x and T are also given. 
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FIG. 8: (Color online) (a) Plots of the probability distribu- 
tion P q versus q/qi for x = 0.2 and T/x = 0.4. o, □, x are 
for L = 10, 8 and 6, respectively. The thick dashed line is 
for the Gaussian distribution that ensues for a paramagnet in 
the macroscopic limit, (b) Same as in (a) but for T = 0.16. 
(c) Same as in (a) but for T — 0.12. Error bars are shown 
wherever they are larger than symbol sizes. 

Ir = I n macroscopic paramagncts, q r is ex- 

pected to be normally distributed, as follows from the 
law of large numbers and the fact that spin-spin cor- 



relation lengths are then finite. On the other hand, 
P q = [5(q r — 1) + 5(q r — l)]/2, where S is the Dirac delta 
function, in a SG phase, according to the droplet picture 
of SGs. 37 Plots of P q vs q r are shown for x = 0.2 in Figs. 
Ek, [8]b, andHt- Clearly, P q (q r ) drifts with system size in 
Fig. [5^,, for T = 0.28. Our results are consistent with 
Pq(lr) — ► (1/tt) ex p( — Qr/ W ) as N ~ * °°; which is in ac- 
cordance with a paramagnetic phase. On the other hand, 
we find for lower temperatures double peaked distribu- 
tions in Figs. [5J} and Fig. [Sfc that are fairly broad and, 
within errors, do not change with N. This is contrary to 
the prediction of the droplet-model theory of SGs. From 
these graphs we conclude that 0.16 < T sg < 0.26 for 
x = 0.2. Analogous plots for x = 0.35 (not shown) give 

0. 30 < T sg < 0.45. 

Results for the scale free quantity A 2 follow. Recall 
that, as explained for A^, A^ — > 7r/2 — 1 a,s N ^ 00 
in the paramagnetic phase, vanishes when there is strong 
long-range order, and goes, at the critical temperature, 
to some intermediate value that is size independent. This 
is as shown in Fig. \Eb for x = 0.7 where curves for various 
values of A" cross at T^p. Figures [5^ and [5k look rather 
similar, because q and m are not qualitatively different 
in the AF phase. This is not so for x < x c , where there 
is no AF order. Figures[5}3 and[5j: for x = 0.6 show that, 
within errors, curves of A^ vs T for different system sizes 
merge (not cross) near T = 0.65, while A^ n increases 
with A^ for all temperatures. Similarly, A^ vs T curves 
merge, for x = 0.65, near T = 0.75 (not shown). Plots 
of A^ vs T/x are shown in Figs. [Hk, Wh>- and[5J: for lower 
concentrations. 

We notice that curves in Figs. [9ji, [9J3, and [9h differ 
only slightly. This follows from the argument given in 
Sec. I, which shows that all physical quantities for three 
dimensional dipolar systems can only be functions of T/x 
for x -C 1. The data points in Fig. [9] show that A^ — >• 
7r/2 — 1 as A^ — s- 00, for T/x > 1, as expected for the 
paramagnetic phase. 

Curves for A^ vs T seem to merge at a lower tempera- 
ture, near T/x = 0.9. However, closer scrutiny shows 
that these curves actually cross, albeit at very small 
glancing angles. This can be appreciated in Figs. [9ji, 
Ek, and[9f, where plots of the ratios A 2 q (L) / A 2 q (A) vs. T 
are given for various values of L, for x = 0.5, x = 0.35, 
and x = 0.2, respectively. Note that the weak depen- 
dence of A^ with system size at low temperatures is in 
accordance with our result that P q (q r ) does not change 
appreciably with system size below T sg . This point is 
further elaborated in Sec. IV Bl 

Following the lead of Refs. and lU who have 
found that £ L /L (defined in Sec. IITC]) crosses at T sq 
and spreads out as T decreases below T sg for the EA 
model in 3D, we next examine how £l/L behaves for 
the PAD model. As pointed out in Sec. I and Table 

1, this has already been done for the PAD model on a 
LiHo ;c Y 1 _ a .Y4 lattice by Kam and Gingras^ As we also 
point out in Sec. I, we aim to explore the behavior of the 
PAD model, not only near T sg , but also deep into the SG 
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FIG. 9: (Color online) (a) Plots of A 2 vs T/x, for x = 0.5. o, 
□ , o, and x are for L = 10, 8, 6 and 4, respectively, (b) Same 
as in (a) but for x = 0.35. (c) Same as in (a) but for x = 0.2. 
(d) Plots of A 2 ,/ A 2 , (4) vs T for x = 0.5. Symbols are as in 
(a), (e) Same as in (d) but for x — 0.35. (f) Same as in (d) 
but for x = 0.2. In panels (a), (b), and (c), all error bars are 
smaller than symbol sizes. 



phase. Recall that £l becomes a true correlation length 
when £l/L <C 1. Then, in the paramagnetic phase, 
£l/L ~ 0(1/ L), therefore decreasing as L increases. At 
T = T sg , ^l/L must become size independent, as ex- 
pected for a scale free quantity. The inferences one can 
make about the nature of the condensed phase from the 
behavior of £l where T < T sg is the subject of Sec. IV Bl 
Without further comment, we next report our results. 
Plots of £l/L versus T/x are shown in Figs. [TOa and 
[TOb for x = 0.35 and 0.2, respectively. Note that curves 
spread out above and below T/x ~ 1. For x = 0.35, 
curves for all L cross at T sg /x = 0.95(5). On the other 
hand, the temperatures where pairs of curves for lengths 
L2 and L\ cross for x = 0.2 decrease as lengths L2 and 



L\ increase (see Fig. [TUb). pointing to a T sg /x < 1.1. 
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FIG. 10: (Color online) (a) Semilog plots of (a) versus 
T/x for x = 0.35, and L = 10 (■), L = 8 (•), L = 6 (♦), and 

L = 4 (a). Dashed line follows from 1/L — > straight line 
extrapolations in the plots shown in Fig. 112b for T < T sg . 
Continuous lines are guides to the eye. (b) Same as in (a) but 
for x = 0.2. All error bars are smaller than symbol sizes. 
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FIG. 11: (Color online) (a) Plots of distributions P q versus 
q/qi for x = 0.2 and the shown values of L and T. Error bars 
are shown only where they are larger than symbol sizes, (b) 
Plots of Tg/x versus JV for the shown values of x. The thick 
dashed line stands for the TV -1 / 2 behavior obtained in Ref. 

M 



V. EXISTENCE AND NATURE OF THE SG 
PHASE 



In this section we examine the numerical results given 
in the previous section. We (i) arrive at values for T sg 
as a function of x, (ii) show that weak long-range order 



is consistent with our results for the SG phase, and (iii) 
draw values for the critical exponent v for various values 
of x. 
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The value of T s , 



Recall first that A 2 vs T curves for different values of 
L are supposed to come together as T approaches T sg 
from above. This behavior is exhibited in Figs. [9K-c. A 
closer view of how such curves actually meet at T = T sg 
is offered in Figs. [HJi-f, where plots of y(L, 4) versus T/x, 
where y(L,L') = A 2 (L)/ A 2 (L'), are shown. One aims 
to find the L — > oo and V — > oo limit of y(L,L') = 1, 
which gives the value of T sg . We find that y(L,L') = 1 
at values of T/x that increase with L and L 1 , which is 
reassuring, because it shows that T sg does not vanish. 
Furthermore, we draw the following lower bounds from 
the plots in Figs. Hf-f, T sg /x > 0.95, 0.8, 0.95, for 
x = 0.5, 0.35, 0.20, respectively. 

We obtain a complementary determination of T sg from 
the intersection of vs T curves. This is as is some- 
times done for the EA4^~— and models. We ob- 
tain, from Fig. [TUJi, T sg /x ~ 0.95 for x = 0.35. In 
Fig. [TUb . we see that ^l/Lvs T curves meet at decreas- 
ingly smaller values of T as L increases. We thus obtain 
Tsg/x < 1.1 for x = 0.2. 

From these two complementary determinations, we ar- 
rive at: T sg /x = 1.0(1) for x < 0.5. 

An aside follows about the result by Snider and Yu? 
that T sg = for x = 0.045, 0.12 or 0.2. This is, of course, 
in clear contradiction with our results. Their conclusions 
come from their work with the Wang-Landau^ variation 
of the MC algorithm. Their evidence is from plots of 
T g versus N, where T g is the temperature at which P q 
becomes flattest. This procedure makes sense because 



.33 



T, 



T ag as N 



-1/2 



oo. They found T g to vanish as N~ 
for several x values, including x = 0.2. We now repeat 
this procedure using our own data, including the ones for 
x = 0.2. In Fig. [TTa we plot the flattest distributions we 
found for x = 0.2 and L = 4, 8, and 10. Note in passing 
that all scaled distributions coincide and have therefore 
the same value of A 2 . Plots of the values of T g /x we 
have obtained for x = 0.5, 0.35, and 0.2 are shown in 
Fig. [Tib . Our data points are in clear contrast to the 
T g ~ A" 1 / 2 trend of Ref. HI, and point to T sg /x ~ 1. 
Whether this disagreement comes from using a different 
Monte Carlo method, or from the unusual definition of q 
in Ref. we do not know. 



B. Marginal behavior 

Here we discuss how various pieces of evidence (in- 
cluding crossings of £,l/L vs T curves) lead us to the 
conclusion that the SG phase of the PAD model behaves 
marginally. That is to say, that (q 2 ) — > and \sg — > oo 
in the macroscopic limit. 

The variation of (q 2 ) with L for various temperatures, 
exhibited in Figs. [7K-c, has already been considered in 
Sec. [TV] For all x < x c , T < T sg , and all system sizes we 
have studied, we find no deviation from (q 2 ) ~ L~( 1+v \ 
Nor do we find any size dependence in P q (q r ). This is il- 
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FIG. 12: (Color online) (a) Semilog plots of £l/L versus 1/L 
for x = 0.35, and T/x = 0.143 (•), T/x = 0.286 (■), T/x = 
0.571 (♦), T/x = 0.857 (T), T/x = 1.00(a), T/x = 1.14 
(0), T/x = 1.43 (□), T/x = 1.71 (o), and T/x = 2.00 (V). 
(b) Same as in (a), but for for x = 0.20, and T/x = 0.300 
(•), T/x = 0.500 (■), T/x = 0.700 (♦), T/x = 1.00 (▼), 
T/x = 1.10(A), T/x = 1.30 (0), T/x = 1.50 (□), T/x = 2.00 
(o), and T/x = 2.50 (V). All errors are: between 2% and 3% 
in (a), and between 2% and 4% in (b), and are thus hidden 
behind the icons. In both (a) and (b), the straight-dashed 
lines give Xr < 1 fitting values, except for T/x — 1.0 in (b), 



for which \r 



3.3 



lustrated in Figs. [5Jd and c, and is in accordance with the 
behavior of the distribution of the magnetization that is 
observecU^ in the condensed phase of the 2D XY model. 
Note that the variation of A 2 with system size is a mea- 
sure of the variation of P q (q r ). The very small changes 
we have observed in A 2 as L varies in the PAD model for 
all T < T sg turn out to be smaller than the correspond- 
ing changes in the XY models This is, of course, in 
marked contrast with the behavior one expects of the cor- 
responding quantity for a strongly ordered system, such 
as the droplet model of SGs or an ordinary fcrromagnct, 
in which A 2 — > in the macroscopic limit of the ordered 
phase. Neither do our results fit into a RSB scenario^ 
in which qi does not vanish as L — > oo and would have 
P q {q r ) changing with system size, since P q {q) is wide and 
does not change with system size in the SG phase. 

We now analyze the data we have for First, we 
outline how we expect £l / L to spread out as T decreases 
below T sg in various SG scenarios. 

(i) Condensate with short range order fluctuations. In 
such a SG phase, (?2 7^ and {4>o4> r ) ~ {4>o){ ( t ) r) would be 
short ranged. This would fit into the droplet model of 
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T/T 



FIG. 13: (Color online) Semilog plots of vs T/x from 

Eq. (|T5J| for the shown values of L. In Eq. (TTZjl . we let 
,4 = 0.67, and 77 = -1 + {T/T sg ) 2 . 



spin glasses^ It then follows straightforwardly from its 
definition [Eq. ©] that £, 2 L /L 2 ~ L d . Here, d = 3, and 
there is nothing in the plots of £l/L vs 1/L, which are 
shown in Figs. [T2k andfl"2b. to suggest that £|/L 2 ~ L 3 
at any nonzero temperature. 

(ii) Condensate with long range order fluctuations. 
Let (A) q be the thermal average of A over all states 
with a given q value. Clearly, (A) = J(A) q P q dq. As- 



sume q2 7^ 0, and J[(<po4> r )q — <l 2 ]Pqdq = G(r), where, 



G(r) 



A 



r d-2+ri ■ 



for r a, where A is a constant. This behavior fits in 
with the RSB picture.— Then, it follows from its defini- 
tion [Eq. ©] that (Hi? ~ L 1+7 >. Recall, from Sec. 
JV\ that 77 ~ -1 + (T/T sg ) 2 in the SG phase. Evidence 
for ££/L 2 ~ L 1+n appears neither in Fig. [P2a nor in 
Fig. Eb. 

(hi) Marginal behavior. Then, (72 = and (<po<Pr) = 
G(r). This is as in the KT theory^ of the 2D XY model. 
It then follows straightforwardly from the definition of 
£l/L that £l/L becomes independent of L for very large 
L. This is precisely the outcome from 1/L — > extrap- 
olations of the straight lines shown in Fig. [T2"k and [T2b 
for all T/x < 1. 

Note also in Figs. 12a and 12b that curves for £l/L 
vs 1/L become steeper as T decreases below T/x ~ 1. 
Now, recall from above that (72 7^ implies £ 2 /L 2 ~ L 
and £,l/L 2 ~ L 1+r > . for short- and long-range fluctua- 
tions from the condensate. Note further that | 1 + 77 | 
decreases as T decreases. This would lead to £l/L vs 
1/L curves which do not become steeper as T decreases 
below T/x ~ 1, which is in clear contradiction with the 
observed behavior. This is an additional piece of evidence 
for quasi long-range order. 

Thus, the most straightforward interpretation of the 
data shown in Figs. \Vlk and!12b leads us to suspect that 




FIG. 14: (Color online) (a) Semilog plots of £,l/L versus 
(T/T sg - 1)L 1/U for x = 0.35, T sg = 0.345, v = 0.95, and 
the shown values of L. (b) Same as in (a) but for x = 0.20, 
T ag = 0.21, v = 0.95, and the shown values of L. Recall that 
(12) scaling is expected only for T/T sg — 1 > 0. In both panels, all 



error bars are somewhat smaller than the icon sizes. 



the SG phase in the PAD model behaves marginally. This 
might seem to be in contradiction to the fact that S^l/L 
curves do cross, as shown in Fig. 1101 and that, as pointed 
out in Ref. lU £l/L vs T curves merge, not cross, for 
the 2D XY model, as T — > T sg from above. (Indeed, 
no crossings occur for even much smaller 2D XY systems 
than the ones for which data points are shown in Ref. l4l| ) . 
We next give a specific example in order to illustrate how 
both merging and spreading out as T decreases below T sg 
can take place, depending on the some details in G(r). 



We first calculate £,l/L from 
I2l) for all r except that G(r) 



)0 r ) = G(r) and Eq. 
1 for all r < 1. To 
proceed, we let A = 0.67 for T < T sg but not too close to 
T = 0, where one expects A = 1. We arc not interested 
here in the T > T sg range, but we nevertheless then let 
A -> Ae~ r ^^, Coo = 7(T/T sg - 1)-", and v = 1, which 
is roughly the value we obtain below (see Sect. IV C|) . We 
make use of 77 = —1 + (T/T sg ) 2 , which we have found in 
Sec. IIV1 Finally, in order to be able to make comparisons 
with our MC results, which we have obtained for periodic 
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boundary conditions, we let in Eq. (|12l) . 

3 

r^Q-'^shi 2 ^)] 1 / 2 , (13) 

a=l 

where Q = ir/L and r = {r\,r2,Ts}- Straightforward nu- 
merical implementation of Eq. ([9]) yields the data points 
that are plotted in Fig. [T2J Note the resemblance be- 
tween Fig. [T3] and Figs. [TOk and 1 10b which follow from 
our MC calculations. 

Merging of £l/L curves at T = T sg as T decreases 
is obtained for all L > 4 if, instead of A = 0.667, we let 
3A = 3-(T/T sg ) 2 . Note that A(T sg ) = 0.667 and A(0) = 
1. If, on the other hand, one lets 3 A = 3 - {T/T sg ) s 
and < s < 0.2, which satisfies the same end point 
conditions, one obtains plots for £l/L vs T which look 
much like the ones shown in Figs. 1101 

To summarize, all our data (including spreading out 
of £,l/L curves as T decreases below T sg ) are consistent 
with marginal behavior in which the correlation length 
diverges at T sg as in a conventional phase transition, but 
weak- long-range order occurs below T sg , as in the 2D XY 
model. 



C. The v exponent 

In accordance with the above, we look for the values of 
v and T sg which best collapse £ t /L vs (T/T sg - l)L x / v 
plots for various values of L into a single curve for tem- 
peratures above T sg . The best results, exhibited in Figs. 
[14k and [Mb . for x = 0.35 and x = 0.20, are obtained 
with T sg /x = 1-0(1) and v = 0.95. Note the data points 
scatter below T sg . This is as expected, and is consistent 
with quasi-long range order in the SG phase, since £l/L 
becomes independent of L then for sufficiently large L. 
Note that, as in the EA model^ 2 - L = 4 seems to be too 
small to scale properly. 



VI. DISCUSSION 

By tempered Monte Carlo calculations, we have stud- 
ied an Ising model on a simple cubic lattice. There are 
only dipole-dipolc interactions. Spins (randomly) occupy 
only a fraction x of all lattice sites. We have calculated 
the entire phase diagram of the system. It is shown in 
Fig. [2] We have also provided strong evidence for the 
existence a SG phase for < x < x c , where x c — 0.65(5). 
The SG transition temperature is given by T sg (x) ~ x. 
We have argued in Sec. I that this result carries over 
into other lattices if (i) x <C 1, and (ii) we replace the 
latter expression for T sg by ksT sg = n^Sd ( sec Table 
I) . How we have arrived a this conclusion is described in 
Sec. IVAl 

We have not dwelt on the applicability of our MC re- 
sults to experiments. That is beyond the scope of this 



paper. We nevertheless make a few comments. Recall 
first that, as we argue in Sec. I, lattice structure is of 
no consequence for very dilute PAD models. Then, T sg 
as well as the temperature T m where the specific heat 
takes its maximum value can only depend (as in the 
MC simulations of Ref. l3lT ) on n^ed (see Table I). We 
notice in Table I values for T sg do not fully comply 
with this rule. In addition, in very dilute LiHo x Yi_j;F4 
systems, T m hardly changes with x*£ There are sev- 
eral sources for the discrepancies between experiments 
on very dilute LiHo a; Yi_ a ;F4 and the PAD model. Quan- 
tum effects seem to play a role in experiments on very di- 
lute LiH0j.Y1_j.F4 systems^ This is not too surprising, 
since tunneling can become relevant when barrier ener- 
gies become overwhelmingly large. However, we do not 
expect small perturbations that bring about tunneling 
and concomitant time dependent effects to have a signifi- 
cant effect on equilibrium properties, which is the subject 
of this paper. In addition, exchange couplings among 
nearest neighbor spin o 31 ' 48 are disregarded in the PAD 
model we study. Note, however that the effect of nearest 
neighbor interactions must vanish as 1 — > 0. Cluster- 
ing of the spatial distribution of dipoles can also lead to 
discrepancies^ None of the above can however account 
for (i) the numerical differences between the MC results 
(see Table I ) of Tarn and Gingras^ 2 - and ours, nor can 
they account for the more serious discrepancy with (ii) 
Ref. HH, which we discuss in some detail in Sec. IVAl 
Numerical (not too large) discrepancies notwithstanding, 
our results support the ones from Tarn and GingraSf 3 ^ 
that the dilute PAD model does have a SG phase. On 
the other hand, for the roots of the discrepancies with ex- 
perimental results (see Table I ) on dilute LiHo a; Y 1 _ :I .F4 
systems, we have no clear picture. 

As for the nature of the SG phase, all of our results 
are consistent with quasi-long-range order. Full details 
are given in Sec. IV Bl We know of no previous study of 
the nature of the SG phase of the PAD model with which 
to compare our results. (Only the critical behavior of a 
PAD model is examined in Ref. [U) On the other hand, 
our conclusion for the PAD model can be compared with 
and one drawn for the EA model in Rcfs. l4£H42l They 
are both based on the behavior of ^/i vs T curves for 
various values of L. The conclusions differ, not so much 
because of the data, but because we have looked at the 
data differently (see Sec. EH and Refs. S^-El). 
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Appendix A: WHY WE DO NOT DO EWALD 
SUMS 



We consider site-diluted systems of Ising magnetic 
dipoles in a cubic box of L 3 sites on a SC lattice. All 
dipoles point along the z axis of the lattice. Each site is 
occupied with probability x. We assume thermal equi- 
librium. We show two things in this appendix. We first 
show that the contribution Ah to the magnetic field h at 
the center of such box, coming from a periodic arrange- 
ment of replicas that span all space beyond the system 
of interest (the "outer space" ) within an arbitrarily large 
cube which is centered on the system of interest, vanishes 
as L — > oo if the system is not in a ferromagnetic phase 
or close to its Curie temperature. More specifically, we 
show that if (si-Sj) — (si){sj) is short ranged, and the 
system is homogeneous (including antiferromagnetically 
ordered states), then 



(A/* 2 







(Al) 



as L — > oo, where (. . .) stands for an average over 
both a canonical ensemble and (site occupation) disor- 
der. Note that we are not imposing the condition that 
(siSj) 2 — (si) 2 (sj) 2 be short ranged, and recall (1) that 
in general J2j( s i s j) ~ ( s i)( s j) = T Xm, where Xm is the 
magnetic susceptibility per site, and (2) that Txm ^$ 1 
for spin glasses. Equation (|Alj) clearly indicates that 
thermodynamic limits can be obtained from Monte Carlo 
calculations for systems of various sizes in which contri- 
butions from the outer space are disregarded. Finally, 
explicit numerical evidence, Fig. [T5J to this effect is also 
given. 



To begin, let h = J2j T ijS] ( Ah = J2j T ijSj) be the 
sum is over all occupied sites within (outside) a cubic 
box of L x L x L sites, centered on i. Therefore, 



Ah 2 



Tin Tim $n Sm 
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FIG. 15: (Color online) Semilog plots of A^l/L vs T/x, where 
the A^l/L is the difference between correlation lengths we 
report in this paper and correlation lengths that obtain when 
Ewald sums are included for x — 0.35, (♦) L — 4 and (•) 
L = 8. These data points follow from averages over 10 4 and 
5 x 10 3 systems samples, for L — 4 and L — 6, respectively. 
The same sample realizations were used for the calculations 
with and without Ewald sums. 



where the sum is over all occupied sites within our system 
of interest. We now replace s n by (s n )+Ss n , and similarly 
for s m , in the equation above. Now, it can be checked 
straightforwardly (i) that £) m /(r m )(s„) = if (s n ) is ei- 
ther independent of n (which would not hold for a ferro- 
magnet with domains) and (ii) that J2 m f( Vm )( Sn ) ~~ * ® 
as L — > oo if (s n ) follows an antiferromagnetic order 
(which, for up and down spins with dipolar interactions 
on a SC lattice, is a checkerboard-like arrangement of up 
and down ferromagnetic columns). Performing thermal 
and disorder averages over the above equation, one then 
obtains, 



(Ah 2 



f(r n )f(r m )(Ss n Ss ri 



(A5) 



as L — > oo. Now, /(r) varies smoothly within the system, 
whence 



(A^)^^[/(r„)] 2 ^(^„^; 



(A6) 



where the double sum is over all occupied sites in the 
outer space. Let 



/(r„) 



E 

j 



£ a a" 



3(z r 



Zj) : 



(A3) 



where Rj is the position of the outer jth box, r n is the 
n-th site's position with respect to the center of the box, 
and the sum is over all outer boxes. Equation (|A2p then 
becomes, 



Ah 2 = y^J(r n )f(r m )s n s T 



(A4) 



if (5s n Ss m ) ~ unless | r„ — r m |<c L. Finally 
]rj/(r„)] 2 = xbe 2 a /L 3 , where b ~ 7.6 if L > 1, as 
follows straightforwardly by numerical integration. Re- 
placement of J2 m ( Ss n$s m ) by Tx m gives Eq. dSU if 
Tx-m is finite. For all the parameters used in our MC 
calculations, we have found that Txm ;$ 1- 

The difference A£l /L between the correlation lengths 
we report and the ones obtained when Ewald sums^ are 
included, for two system sizes, are exhibited in Fig. [THl 
The same sample realizations were used for the calcula- 
tions with and without Ewald sums. This explains why 
we can show in Fig. I15l valucs for A£l/L that are smaller 
than the statistical errors given for £^ /L (see Fig. [T2^) for 
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L = 6. The results are clearly consistent with a A£l/L that vanishes in the thermodynamic limit. 
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